#Project: TV polarization
#Author: Josh Ryan and Jeff Lyons
#Working Date: Winter 23
#Task: Use Imai, Kim, Wang PanelMatch analysis--find best matching algorithm
#Related files: leg_level_forstata_v6_102524.dta is loaded


library(foreign)
library(effects)
library(sciplot)
library(sandwich)
library(lmtest)   
library(BMS)
library(rJava)
library(lattice)
library(plyr)
library(lme4)
library(gmodels)
library(MASS)
library(stringr)
library("XLConnect")
library(reshape2)
library(xlsx)
library(dplyr)
library(doBy)
library(Hmisc)
library(stargazer)
library(xtable)
library(dplyr)
library(xtable)
library(foreign)
library(dotwhisker)#https://cran.r-project.org/web/packages/dotwhisker/vignettes/dotwhisker-vignette.html
library(berryFunctions) #for the insert row function
library(ggplot2)
library(gridExtra)
library(wnominate)
library(tidyr)
library(wnominate)
library(oc)
library(haven)
library(PanelMatch)
library(gridGraphics)
library(grid)


####NOTE: THESE BALANCE TESTS WERE ORIGINALLY CONDUCTED USING PANELMATCH 1.0.0. THE NEWEST VERSION OF THE PACKAGE USES DIFFERENT CODE FROM THAT HERE.

dat<-read_dta("leg_level_forstata_v6_102524.dta")


###Figure I2-I4

#Have to keep making as dataframe because of the way Panelmatch works
dat<-as.data.frame(dat)

#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.np<-dat[c("np_score", "treatment", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "distid1")]

#have to constantly be making these as integers, super annoying
dat.np<-as.data.frame(dat.np)
dat.np$distid1<-as.integer(dat.np$distid1)
dat.np$year<-as.integer(dat.np$year)

dat.np0<-subset(dat.np, party_id==0)
dat.np1<-subset(dat.np, party_id==1)

#check balance improvement first run matching results, then look at balance improvement
######np score Dems

##No refinement method
results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                         treatment = "treatment", refinement.method = "none", 
                         data = dat.np0, match.missing = TRUE, 
                         covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                                  I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                         + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                           I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                         size.match = 10, qoi = "att" ,outcome.var = "np_score",
                         lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.none, data = dat.np0)
summary(PE.results)
plot(PE.results)

##Mahalnobis Distance
results.mah <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                          treatment = "treatment", refinement.method = "mahalanobis", 
                          data = dat.np0, match.missing = TRUE, 
                          covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                            I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                          + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                            I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                          size.match = 10, qoi = "att" ,outcome.var = "np_score",
                          lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.mah, data = dat.np0)
summary(PE.results)
plot(PE.results)

#Testing different balance
#####
##Proposenity Score Matching
results.ps <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                          treatment = "treatment", refinement.method = "ps.match", 
                          data = dat.np0, match.missing = TRUE, 
                          covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                            I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                          + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                            I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                          size.match = 10, qoi = "att" ,outcome.var = "np_score",
                          lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.np0)
summary(PE.results)
plot(PE.results)

#Propensity Score Weighting
results.weight <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                          treatment = "treatment", refinement.method = "ps.weight", 
                          data = dat.np0, match.missing = TRUE, 
                          covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                            I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                          + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                            I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                          size.match = 10, qoi = "att" ,outcome.var = "np_score",
                          lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.weight, data = dat.np0)
summary(PE.results)
plot(PE.results)

#covariate balance propsensity score matching
results.cbps <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                             treatment = "treatment", refinement.method = "CBPS.match", 
                             data = dat.np0, match.missing = TRUE, 
                             covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                               I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                             + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                               I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                             size.match = 10, qoi = "att" ,outcome.var = "np_score",
                             lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.cbps, data = dat.np0)
summary(PE.results)
plot(PE.results)

results.cbpsweight <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                           treatment = "treatment", refinement.method = "CBPS.weight", 
                           data = dat.np0, match.missing = TRUE, 
                           covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "np_score",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.cbpsweight, data = dat.np0)
summary(PE.results)
plot(PE.results)

#doesn't converge
# results.psmsmweight <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
#                                  treatment = "treatment", refinement.method = "ps.msm.weight", 
#                                  data = dat.np0, match.missing = TRUE, 
#                                  covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
#                                    I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
#                                  + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
#                                    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
#                                  size.match = 10, qoi = "att" ,outcome.var = "np_score",
#                                  lead = 0:4, forbid.treatment.reversal = TRUE)
# 
# PE.results <- PanelEstimate(sets = results.psmsmweight, data = dat.np0)
# summary(PE.results)
# plot(PE.results)


# results.cbpspsmsmweight <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
#                                   treatment = "treatment", refinement.method = "CBPS.msm.weight", 
#                                   data = dat.np0, match.missing = TRUE, 
#                                   covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
#                                     I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
#                                   + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
#                                     I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
#                                   size.match = 10, qoi = "att" ,outcome.var = "np_score",
#                                   lead = 0:4, forbid.treatment.reversal = TRUE)
# 
# PE.results <- PanelEstimate(sets = results.cbpspsmsmweight, data = dat.np0)
# summary(PE.results)
# plot(PE.results)



#A dot below the 45 degree line implies that the standardized mean balance is improved after the refinement for a particular time-varying covariate.
 #Testing different matching techniques
########################
par(mfrow=c(2,3))

mah.balance<-balance_scatter(matched_set_list = list(results.mah$att),
                main = "Mahalanobis Distance",
                data = dat.np0,
                covariates = c("np_score", "majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))

ps.balance<-balance_scatter(matched_set_list = list(results.ps$att),
                                       main = "Propensity Score",
                data = dat.np0,
                covariates = c("np_score", "majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))

psweight.balance<-balance_scatter(matched_set_list = list(results.weight$att),
                                main = "Propensity Score Weighting",
                data = dat.np0,
                covariates = c("np_score", "majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))

cbps.balance<-balance_scatter(matched_set_list = list(results.cbps$att),
                                    main = "Covariate Balance Prop. Score Matching",
                                    data = dat.np0,
                                    covariates = c("np_score", "majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))


cbpsweight.balance<-balance_scatter(matched_set_list = list(results.cbpsweight$att),
                                main = "Covariate Balance Propensity Score Weighting",
                                data = dat.np0,
                                covariates = c("np_score", "majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))


get_covariate_balance(results.mah$att, dat.np0, covariates = c("majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                            "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.ps$att, dat.np0, covariates = c("majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                          "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.weight$att, dat.np0, covariates = c("majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                           "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.cbps$att, dat.np0, covariates = c("majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                  "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.cbpsweight$att, dat.np0, covariates = c("majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                 "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

###############mahalanobis performs the best



######np score GOP

##No refinement method
results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                           treatment = "treatment", refinement.method = "none", 
                           data = dat.np1, match.missing = TRUE, 
                           covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "np_score",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.none, data = dat.np1)
summary(PE.results)
plot(PE.results)

##Mahalnobis Distance
results.mah <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                          treatment = "treatment", refinement.method = "mahalanobis", 
                          data = dat.np1, match.missing = TRUE, 
                          covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                            I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                          + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                            I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                          size.match = 10, qoi = "att" ,outcome.var = "np_score",
                          lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.mah, data = dat.np1)
summary(PE.results)
plot(PE.results)

##Proposenity Score Matching
results.ps <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                         treatment = "treatment", refinement.method = "ps.match", 
                         data = dat.np1, match.missing = TRUE, 
                         covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                           I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                         + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                           I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                         size.match = 10, qoi = "att" ,outcome.var = "np_score",
                         lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.np1)
summary(PE.results)
plot(PE.results)

#Propensity Score Weighting
results.weight <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                             treatment = "treatment", refinement.method = "ps.weight", 
                             data = dat.np1, match.missing = TRUE, 
                             covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                               I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                             + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                               I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                             size.match = 10, qoi = "att" ,outcome.var = "np_score",
                             lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.weight, data = dat.np1)
summary(PE.results)
plot(PE.results)

#covariate balance propsensity score matching
results.cbps <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                           treatment = "treatment", refinement.method = "CBPS.match", 
                           data = dat.np1, match.missing = TRUE, 
                           covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "np_score",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.cbps, data = dat.np1)
summary(PE.results)
plot(PE.results)

results.cbpsweight <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                                 treatment = "treatment", refinement.method = "CBPS.weight", 
                                 data = dat.np1, match.missing = TRUE, 
                                 covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                                   I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                                 + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                                   I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                                 size.match = 10, qoi = "att" ,outcome.var = "np_score",
                                 lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.cbpsweight, data = dat.np1)
summary(PE.results)
plot(PE.results)


#A dot below the 45 degree line implies that the standardized mean balance is improved after the refinement for a particular time-varying covariate.

########################
par(mfrow=c(2,3))

mah.balance<-balance_scatter(matched_set_list = list(results.mah$att),
                             main = "Mahalanobis Distance",
                             data = dat.np1,
                             covariates = c("np_score", "majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))

ps.balance<-balance_scatter(matched_set_list = list(results.ps$att),
                            main = "Propensity Score",
                            data = dat.np1,
                            covariates = c("np_score", "majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))

psweight.balance<-balance_scatter(matched_set_list = list(results.weight$att),
                                  main = "Propensity Score Weighting",
                                  data = dat.np1,
                                  covariates = c("np_score", "majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))

cbps.balance<-balance_scatter(matched_set_list = list(results.cbps$att),
                              main = "Covariate Balance Prop. Score Matching",
                              data = dat.np1,
                              covariates = c("np_score", "majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))


cbpsweight.balance<-balance_scatter(matched_set_list = list(results.cbpsweight$att),
                                    main = "Covariate Balance Propensity Score Weighting",
                                    data = dat.np1,
                                    covariates = c("np_score", "majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))



get_covariate_balance(results.mah$att, dat.np1, covariates = c("majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                               "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.ps$att, dat.np1, covariates = c("majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                              "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.weight$att, dat.np1, covariates = c("majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                  "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.cbps$att, dat.np1, covariates = c("majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.cbpsweight$att, dat.np1, covariates = c("majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                      "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

#cbsp is best


#########################
#SLES



#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.sles<-dat[c("sles", "treatment", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "distid1")]

dat.sles$distid1<-as.integer(dat.np$distid1)
dat.sles$year<-as.integer(dat.np$year)
dat.sles$sles<-as.integer(dat.sles$sles)


#check balance improvement first run matching results, then look at balance improvement

##No refinement method
results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                           treatment = "treatment", refinement.method = "none", 
                           data = dat.sles, match.missing = TRUE, 
                           covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "sles",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.none, data = dat.sles)
summary(PE.results)
plot(PE.results)

##Mahalnobis Distance
results.mah <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                          treatment = "treatment", refinement.method = "mahalanobis", 
                          data = dat.sles, match.missing = TRUE, 
                          covs.formula = ~ I(lag(majority, 1:4))+I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                            I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                          + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                            I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                          size.match = 10, qoi = "att" ,outcome.var = "sles",
                          lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.mah, data = dat.sles)
summary(PE.results)
plot(PE.results)

##Proposenity Score Matching
results.ps <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                         treatment = "treatment", refinement.method = "ps.match", 
                         data = dat.sles, match.missing = TRUE, 
                         covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4))+ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                           I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                         + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                           I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                         size.match = 10, qoi = "att" ,outcome.var = "sles",
                         lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.sles)
summary(PE.results)
plot(PE.results)

#Propensity Score Weighting
results.weight <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                             treatment = "treatment", refinement.method = "ps.weight", 
                             data = dat.sles, match.missing = TRUE, 
                             covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                               I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                             + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                               I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                             size.match = 10, qoi = "att" ,outcome.var = "sles",
                             lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.weight, data = dat.sles)
summary(PE.results)
plot(PE.results)

#covariate balance propsensity score matching
results.cbps <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                           treatment = "treatment", refinement.method = "CBPS.match", 
                           data = dat.sles, match.missing = TRUE, 
                           covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "sles",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.cbps, data = dat.sles)
summary(PE.results)
plot(PE.results)

results.cbpsweight <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                                 treatment = "treatment", refinement.method = "CBPS.weight", 
                                 data = dat.sles, match.missing = TRUE, 
                                 covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                                   I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                                 + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                                   I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                                 size.match = 10, qoi = "att" ,outcome.var = "sles",
                                 lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.cbpsweight, data = dat.sles)
summary(PE.results)
plot(PE.results)

#doesn't converge
# results.psmsmweight <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
#                                  treatment = "treatment", refinement.method = "ps.msm.weight", 
#                                  data = dat.sles, match.missing = TRUE, 
#                                  covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
#                                    I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
#                                  + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
#                                    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
#                                  size.match = 10, qoi = "att" ,outcome.var = "sles",
#                                  lead = 0:4, forbid.treatment.reversal = TRUE)
# 
# PE.results <- PanelEstimate(sets = results.psmsmweight, data = dat.sles)
# summary(PE.results)
# plot(PE.results)


# results.cbpspsmsmweight <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
#                                   treatment = "treatment", refinement.method = "CBPS.msm.weight", 
#                                   data = dat.sles, match.missing = TRUE, 
#                                   covs.formula = ~ I(lag(majority, 1:4))+I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
#                                     I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
#                                   + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
#                                     I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
#                                   size.match = 10, qoi = "att" ,outcome.var = "sles",
#                                   lead = 0:4, forbid.treatment.reversal = TRUE)
# 
# PE.results <- PanelEstimate(sets = results.cbpspsmsmweight, data = dat.sles)
# summary(PE.results)
# plot(PE.results)



#A dot below the 45 degree line implies that the standardized mean balance is improved after the refinement for a particular time-varying covariate.

########################
par(mfrow=c(2,3))

mah.balance<-balance_scatter(matched_set_list = list(results.mah$att),
                             main = "Mahalanobis Distance",
                             data = dat.sles,
                             covariates = c("sles", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))

ps.balance<-balance_scatter(matched_set_list = list(results.ps$att),
                            main = "Propensity Score",
                            data = dat.sles,
                            covariates = c("sles", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))

psweight.balance<-balance_scatter(matched_set_list = list(results.weight$att),
                                  main = "Propensity Score Weighting",
                                  data = dat.sles,
                                  covariates = c("sles", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))

cbps.balance<-balance_scatter(matched_set_list = list(results.cbps$att),
                              main = "Covariate Balance Prop. Score Matching",
                              data = dat.sles,
                              covariates = c("sles", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))


cbpsweight.balance<-balance_scatter(matched_set_list = list(results.cbpsweight$att),
                                    main = "Covariate Balance Propensity Score Weighting",
                                    data = dat.sles,
                                    covariates = c("sles", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))



get_covariate_balance(results.mah$att, dat.sles, covariates = c("majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                               "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.ps$att, dat.sles, covariates = c("majority", "veto", "party_id", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                              "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.weight$att, dat.sles, covariates = c("majority", "veto", "party_id", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                  "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.cbps$att, dat.sles, covariates = c("majority", "veto", "party_id", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.cbpsweight$att, dat.sles, covariates = c("majority", "veto", "party_id", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                      "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))


#####ps works best


########party loyalty


#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.loyalty<-dat[c("party_loyalty", "treatment", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "distid1")]

dat.loyalty$distid1<-as.integer(dat.np$distid1)
dat.loyalty$year<-as.integer(dat.np$year)
dat.loyalty$sles<-as.integer(dat.loyalty$party_loyalty)


#check balance improvement first run matching results, then look at balance improvement

##No refinement method
results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                           treatment = "treatment", refinement.method = "none", 
                           data = dat.loyalty, match.missing = TRUE, 
                           covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "party_loyalty",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.none, data = dat.loyalty)
summary(PE.results)
plot(PE.results)

##Mahalnobis Distance
results.mah <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                          treatment = "treatment", refinement.method = "mahalanobis", 
                          data = dat.loyalty, match.missing = TRUE, 
                          covs.formula = ~ I(lag(majority, 1:4))+I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                            I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                          + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                            I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                          size.match = 10, qoi = "att" ,outcome.var = "party_loyalty",
                          lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.mah, data = dat.loyalty)
summary(PE.results)
plot(PE.results)

##Proposenity Score Matching
results.ps <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                         treatment = "treatment", refinement.method = "ps.match", 
                         data = dat.loyalty, match.missing = TRUE, 
                         covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4))+ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                           I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                         + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                           I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                         size.match = 10, qoi = "att" ,outcome.var = "party_loyalty",
                         lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.loyalty)
summary(PE.results)
plot(PE.results)

#Propensity Score Weighting
results.weight <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                             treatment = "treatment", refinement.method = "ps.weight", 
                             data = dat.loyalty, match.missing = TRUE, 
                             covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                               I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                             + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                               I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                             size.match = 10, qoi = "att" ,outcome.var = "party_loyalty",
                             lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.weight, data = dat.loyalty)
summary(PE.results)
plot(PE.results)

#covariate balance propsensity score matching
results.cbps <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                           treatment = "treatment", refinement.method = "CBPS.match", 
                           data = dat.loyalty, match.missing = TRUE, 
                           covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "party_loyalty",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.cbps, data = dat.loyalty)
summary(PE.results)
plot(PE.results)

results.cbpsweight <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                                 treatment = "treatment", refinement.method = "CBPS.weight", 
                                 data = dat.loyalty, match.missing = TRUE, 
                                 covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                                   I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                                 + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                                   I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                                 size.match = 10, qoi = "att" ,outcome.var = "party_loyalty",
                                 lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.cbpsweight, data = dat.loyalty)
summary(PE.results)
plot(PE.results)

#doesn't converge
# results.psmsmweight <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
#                                  treatment = "treatment", refinement.method = "ps.msm.weight", 
#                                  data = dat.loyalty, match.missing = TRUE, 
#                                  covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
#                                    I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
#                                  + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
#                                    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
#                                  size.match = 10, qoi = "att" ,outcome.var = "party_loyalty",
#                                  lead = 0:4, forbid.treatment.reversal = TRUE)
# 
# PE.results <- PanelEstimate(sets = results.psmsmweight, data = dat.loyalty)
# summary(PE.results)
# plot(PE.results)


# results.cbpspsmsmweight <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
#                                   treatment = "treatment", refinement.method = "CBPS.msm.weight", 
#                                   data = dat.loyalty, match.missing = TRUE, 
#                                   covs.formula = ~ I(lag(majority, 1:4))+I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
#                                     I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
#                                   + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
#                                     I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
#                                   size.match = 10, qoi = "att" ,outcome.var = "party_loyalty",
#                                   lead = 0:4, forbid.treatment.reversal = TRUE)
# 
# PE.results <- PanelEstimate(sets = results.cbpspsmsmweight, data = dat.loyalty)
# summary(PE.results)
# plot(PE.results)



#A dot below the 45 degree line implies that the standardized mean balance is improved after the refinement for a particular time-varying covariate.

########################
par(mfrow=c(2,3))

mah.balance<-balance_scatter(matched_set_list = list(results.mah$att),
                             main = "Mahalanobis Distance",
                             data = dat.loyalty,
                             covariates = c("party_loyalty", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))

ps.balance<-balance_scatter(matched_set_list = list(results.ps$att),
                            main = "Propensity Score",
                            data = dat.loyalty,
                            covariates = c("party_loyalty", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))

psweight.balance<-balance_scatter(matched_set_list = list(results.weight$att),
                                  main = "Propensity Score Weighting",
                                  data = dat.loyalty,
                                  covariates = c("party_loyalty", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))

cbps.balance<-balance_scatter(matched_set_list = list(results.cbps$att),
                              main = "Covariate Balance Prop. Score Matching",
                              data = dat.loyalty,
                              covariates = c("party_loyalty", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))


cbpsweight.balance<-balance_scatter(matched_set_list = list(results.cbpsweight$att),
                                    main = "Covariate Balance Propensity Score Weighting",
                                    data = dat.loyalty,
                                    covariates = c("party_loyalty", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))



get_covariate_balance(results.mah$att, dat.loyalty, covariates = c("majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.ps$att, dat.loyalty, covariates = c("majority", "veto", "party_id", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                               "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.weight$att, dat.loyalty, covariates = c("majority", "veto", "party_id", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                   "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.cbps$att, dat.loyalty, covariates = c("majority", "veto", "party_id", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                 "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.cbpsweight$att, dat.loyalty, covariates = c("majority", "veto", "party_id", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                       "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

#ps again


###########Figure A2


#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.np<-dat[c("np_score", "treatment", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "distid1")]

#have to constantly be making these as integers, super annoying
dat.np<-as.data.frame(dat.np)
dat.np$distid1<-as.integer(dat.np$distid1)
dat.np$year<-as.integer(dat.np$year)

dat.np0<-subset(dat.np, party_id==0)
dat.np1<-subset(dat.np, party_id==1)

#check balance improvement first run matching results, then look at balance improvement
######np score Dems

##No refinement method
results.none.demideo <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                                   treatment = "treatment", refinement.method = "none", 
                                   data = dat.np0, match.missing = TRUE, 
                                   covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                                     I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                                   + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                                     I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                                   size.match = 10, qoi = "att" ,outcome.var = "np_score",
                                   lead = 0:4, forbid.treatment.reversal = TRUE)

#covariate balance propsensity score matching
results.cbps.demideo <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                                   treatment = "treatment", refinement.method = "CBPS.match", 
                                   data = dat.np0, match.missing = TRUE, 
                                   covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                                     I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                                   + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                                     I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                                   size.match = 10, qoi = "att" ,outcome.var = "np_score",
                                   lead = 0:4, forbid.treatment.reversal = TRUE)

par(mfrow=c(2,2))

cbps.balance<-balance_scatter(matched_set_list = list(results.cbps.demideo$att),
                              main = "Covariate Balance Prop. Score Matching",
                              data = dat.np0,
                              covariates = c("np_score", "majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))


######np score GOP

##No refinement method
results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                           treatment = "treatment", refinement.method = "none", 
                           data = dat.np1, match.missing = TRUE, 
                           covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "np_score",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

results.cbps <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                           treatment = "treatment", refinement.method = "CBPS.match", 
                           data = dat.np1, match.missing = TRUE, 
                           covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "np_score",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

cbps.balance<-balance_scatter(matched_set_list = list(results.cbps$att),
                              main = "Covariate Balance Prop. Score Matching",
                              data = dat.np1,
                              covariates = c("np_score", "majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))




#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.sles<-dat[c("sles","party_id", "treatment", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "distid1")]

dat.sles$distid1<-as.integer(dat.np$distid1)
dat.sles$year<-as.integer(dat.np$year)
dat.sles$sles<-as.integer(dat.sles$sles)


#check balance improvement first run matching results, then look at balance improvement

##No refinement method
results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                           treatment = "treatment", refinement.method = "none", 
                           data = dat.sles, match.missing = TRUE, 
                           covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "sles",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

results.ps <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                         treatment = "treatment", refinement.method = "ps.match", 
                         data = dat.sles, match.missing = TRUE, 
                         covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                           I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                         + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                           I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                         size.match = 10, qoi = "att" ,outcome.var = "sles",
                         lead = 0:4, forbid.treatment.reversal = TRUE)

ps.balance<-balance_scatter(matched_set_list = list(results.ps$att),
                            main = "Propensity Score Matching",
                            data = dat.sles,
                            covariates = c("sles", "majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "party_id"))




#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.loyalty<-dat[c("party_loyalty", "treatment", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "distid1")]

dat.loyalty$distid1<-as.integer(dat.np$distid1)
dat.loyalty$year<-as.integer(dat.np$year)
dat.loyalty$sles<-as.integer(dat.loyalty$party_loyalty)


#check balance improvement first run matching results, then look at balance improvement

##No refinement method
results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                           treatment = "treatment", refinement.method = "none", 
                           data = dat.loyalty, match.missing = TRUE, 
                           covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "party_loyalty",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

results.ps <- PanelMatch(lag = 4, time.id = "year", unit.id = "distid1", 
                         treatment = "treatment", refinement.method = "ps.match", 
                         data = dat.loyalty, match.missing = TRUE, 
                         covs.formula = ~ I(lag(majority, 1:4)) +I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                           I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                         + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                           I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                         size.match = 10, qoi = "att" ,outcome.var = "party_loyalty",
                         lead = 0:4, forbid.treatment.reversal = TRUE)


ps.balance<-balance_scatter(matched_set_list = list(results.ps$att),
                            main = "Propensity Score Matching",
                            data = dat.loyalty,
                            covariates = c("party_loyalty","party_id", "majority", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro"))
